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Abstract. 

We assess the impact of non-thermally shock-accelerated particles on the 
magnetohydrodynamic (MHD) jump conditions of relativistic shocks. The adiabatic 
constant is calculated directly from first principle particle-in-cell simulation data, 
enabling a semi-kinetic approach to improve the standard fluid model and allowing 
for an identification of the key parameters that define the shock structure. We find 
that the evolving upstream parameters have a stronger impact than the corrections 
due to non-thermal particles. We find that the decrease of the upstream bulk speed 
yields deviations from the standard MHD model up to 10%. Furthermore, we obtain 
a quantitative definition of the shock transition region from our analysis. For Weibel- 
mediated shocks the inclusion of a magnetic field in the MHD conservation equations 
is addressed for the first time. 
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1. Introduction 

Shocks are common in the universe and a topic of high interest due to their importance 
in the acceleration of high-energy particles and the subsequent generation of radiation. 
The most prominent examples are the non-relativistic shocks in supernovae, which can 
provide an efficient acceleration of cosmic rays inside our galaxy [l], and relativistic 
shocks in gamma-ray bursts (GRB) |2|. A clear understanding of the shock properties 
and their connection to the structure of the fields and the distribution function of the 
particles is of critical importance to understand and to model many of the scenarios. In 
particular, as laboratory experiments start to explore in detail these conditions [3]-[6] and 
numerical simulations can capture many of the details of these structures (7 16 



more 



detailed theoretical models are also required to explain and to predict the properties of 



relativistic shocks in different contexts 17-19 



The theoretical models to describe the shock properties are based on the 
hydrodynamic jump conditions, and assume a steady state, neglecting the involved 



kinetics. In particular, Blandford and McKee 20 considered strong shocks, which 
appear if either the upstream is cold and the energy per particle stays unchanged or if 
the upstream is ultra-relativistic, so that the rest mass energy can be neglected. In the 
latter case, energy and pressure are connected by the equation of state p = e/3. However, 
and due to the interaction with self-consistent fields in the shock, the particles can be 
trapped and accelerated in the shock, forming the characteristic high-energy tail in the 



distribution function, which has been recently reported in simulations (e.g. 11 12]). 
The standard model of the hydrodynamic jump conditions assumes thermal spectra, 
neglecting the infiuence of accelerated particles. If the non-thermal tail is strong and 
the actual particle distribution deviates from such a spectrum, the pressure and energy 
densities in the downstream vary as well and lead to a modification of the steady state 
conditions, which can be mathematically expressed by a modification of the adiabatic 
constant. 

In this paper we address the effect of such deviations and derive the jump 
conditions based on the actual particle distribution in the shock. In particular, we 
focus on the effects on the shock speed and the density compression ratio which are 
the key parameters for determining the shock dynamics and energy transport. We 
start our analysis with a generalization of the theory for the shock jump conditions 
for an upstream population with non-zero temperature and discuss the impact of 
deviations from the idealized contributing parameters on the jump conditions. The 
theoretical predictions are then compared with fully self-consistent particle-in-cell (PIC) 
simulations. Our analysis demonstrates that the modification of the downstream 
adiabatic constant due to the development of the non-thermal tail as previously reported 
(8|[ll||l2 can have a strong impact, but the decrease of the bulk Lorentz factor directly 
in front of the shock has the dominant infiuence on the jump conditions. Theory and 
simulations can be matched for a well-defined shock transition region, thus contributing 
to identify the different shock regimes. 
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The analysis has been done for a pure electron-positron plasma, as the expected 
effects on the adiabatic constant are qualitatively the same as for electron-ion plasmas 
if the plasma is initially unmagnetized. An initial magnetization suppresses the non- 
thermal acceleration in pure pair plasmas, and the role of ions becomes then important 



in this context 21 



2. Theoretical model 



The starting point for the derivation of the shock jump conditions 20 are the 
conservation equations for mass, momentum and energy. We perform our calculations 
in the downstream rest frame in order to match the configuration of the simulations 
(see next section). In the standard approach, the one-dimensional strong shock 
approximation, the upstream is considered to be cold {pi = 0) and the contributions 
from the self-consistently generated magnetic fields are neglected. Here we include both 



contributions and follow the formalism of 22 , where quantities with a single index Qi 
are measured in their own rest frame and quantities with double indices Qij denote the 
value of species i in the rest frame j. Throughout the paper, indices 1, 2, s refer to 
the upstream, downstream, shock frame, respectively. Thus, the conservation equations 
read 

muis = n2U2s (1) 

PlsBis = P2sB2s (2) 

7is/ii(l + cTi) = 72s/i2(l + 0-2) (3) 

2/3f^ riiuis 2/3|^ ^2^25 

with Uis = ■jisf^is, where ■jis denotes the Lorentz factor, f3is = Vis/c where Vis is the 
bulk velocity, cxj = Bf^/ {Aiinifii'yfJ is the magnetization, where Bis is the transverse 

magnetic field, rij is the plasma density, and /Xj = 1 H — ^— is the specific enthalpy. 

The adiabatic constant Fj is defined by the relation between the energy density Cj and 
pressure density Pi = (Fj — l)(ej — pi) with rest mass density pi = riimc^ . 

We start our analysis by considering the case where the magnetic field contribution 
can be neglected, which is the standard approach for initially unmagnetized shocks 



10 , 20 . We will later discuss the influence of the self-generated magnetic fields on 
the jump conditions in the long time evolution of the shock. The shock speed can be 
determined by performing a Lorentz transformation into the downstream frame and 
combining equations ([T])-(|4]), yielding 

. (F2 - l)(7i2//i - 1) , s 

/^iV7i2 - 1 

which depends only on the upstream Lorentz factor 712, the downstream adiabatic 
constant F2 and the upstream enthalpy pi. A non-zero upstream pressure {pi > 1) 
increases the shock speed. This effect is weaker the higher the upstream Lorentz factor 
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is and approaches the strong shock approximation for yUi = 1 [9j. The density ratio is 
given by 

^ = 1 + ^ = 1+ (7l2 - /gN 

ni2 Ps2 7i2(r2 - l)(7i2Aii - 1) 

and is decreased if the upstream pressure is taken into account. The deviations 
associated with non-thermal tails will have an impact on the adiabatic constant 
In order to assess the influence of small deviations of the adiabatic constant r2 to the 
typically considered adiabatic constant of an ideal gas F^, we rewrite the adiabatic 
constant as r2 = + 6T2 where 6T2 <^ F^. The shock speed is now given by 

A, = + ^2Bpf^sr, (7) 

/^iV7i2 - 1 

and, therefore, the correction of the adiabatic constant increases the shock speed by an 
amount of the order of (5F2 for a highly relativistic upstream. The density ratio 

^ r4_ (7?2 - /g^ 

ni2 ni2 7i2(r2 - l)2(7i2/ii - 1) ^ 
is decreased when the correction of the adiabatic constant is included. Typically, an 
adiabatic constant Fj = 3/2 for 2D and 4/3 for 3D is used to verify the jump conditions 
of relativistic shocks, e.g. |l2j; therefore, the corrections to the density ratio are of the 
order of 4 5F2 in 2D and 9 5F2 in 3D for a relativistic upstream flow. 

Deviations in the Lorentz factor of the flows can also affect the shock jump 
conditions. Following the previous approach, we define the Lorentz factor of the 
upstream flow as 712 = 712 — ^712, where 7^2 is the initial Lorentz factor of the upstream 
flow and 5712 its deviation. An increase of 5712 reduces the shock speed and enhances the 



density ratio according to the Taylor expansion of the jump conditions (see Appendix 
0. 

The effect of the upstream pressure and of deviations of the adiabatic constant and 
upstream Lorentz factor in the density ratio are illustrated in Figure [T| summarizing 
the previous findings and illustrating the stronger impact of the change in the upstream 
Lorentz factor. 

3. Numerical simulations 

In order to address the effect of the different parameters in realistic scenarios, where 
the shock structure evolves in time, in a self-consistent manner, we have performed fully 



relativistic simulations of the shock formation and propagation with OSIRIS 2.0 23 24 
By using a fully kinetic model, the macroscopic quantities describing the shock structure 
can be calculated directly from the kinetic quantities and compared with our theoretical 
model. For a given distribution function /(p) obtained from the simulation data, the 
energy and pressure densities are calculated in the local rest frame of the fluid as 

d'pifip) (9) 



nmc^ 
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Figure 1. Effect of the upstream pressure, and of deviations of the adiabatic constant 
and upstream Lorentz factor in the density ratio. The increase of the upstream pressure 
and the slowdown of the flow increase the density ratio, whereas deviations on the 
adiabatic constant decrease the density ratio. All curves are plotted for 75*2 = 20. 
Black lines correspond to /ii = 1, red lines to /ii = 2, solid lines to r2 = 1.5 and 
dashed lines to r2 = 1.52. 



p:=r^= / d'p^f ip) (10) 



nmc^ J " 7 

with 7 = a/1 + p^. Note that the integrals reduce to double integrals for the 2D case. 
The adiabatic constant can then be calculated from the previously mentioned relation 
between energy and pressure densities as F = 1 +p/ (e — 1). For a relativistic Maxwellian 
/(7) = Cexp(— 7/A7) the adiabatic constant yields F2D = (2 + 3A7)/(1 + 2A7) for 
the 2D case and T^d = I + A7/PA7 - 1 + Ki{A-f-^)/K2iA-f-^)] for a 3D geometry 
with Kn{x) the modified Bessel functions of the second kind. The limiting values are 
F2D = 2 for A7 0, F2D = 3/2 for A7 ^ 00 and Fg^ = 5/3 for A7 ^ 0, Tsd = 4/3 
for A7 — )■ 00. 

We can immediately observe that, assuming a full thermalization of the upstream 
flow in the downstream, with a spread A7 = (712 — l)/2 [lO], the density ratio equation 
(|6| reduces to n2/ni2 = 3 in 2D, independent of the initial upstream Lorentz factor. 
Even for a highly relativistic flow, the deviations arising from the correction of the 
adiabatic constant can be noticeable, for instance n2/ni2 = 3.1 for 712 = 20 and 
f^2l'n\2 = 3.13 for 7i2 = 15 [8|. 

In reality, a more complex distribution function of the particles is expected due to 
the accelerated particle component. Previous results found the best fit for a Maxwellian 
bulk plus a power-law tail 

fil) =^"^^ =Ciexp [-7/A7] 

+ C27""-imin{l,exp[-(7-7e„t)/A7e„t]} (11) 

with C2 = for 7 < 7mm (lO). The contribution of this modified distribution to 
the macroscopic shock properties can be now addressed for the first time, using the 
self-consistent particle distribution from the simulations. The cumbersome analytical 
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expressions for the energy and pressure densities are presented in Appendix B Large 
effects are expected for a very strong tail, which is given by a small 7mm in combination 
with a small a. In the case of relativistic shocks, these parameters are such that the 
contribution from the tail is weak. 

In our simulations, the relativistic shock is created by injecting a charge neutral 
electron-positron beam with an isotropic thermal spread of 10~^c and bulk Lorentz factor 
7°2 = 20 along the negative Xi direction. The particles are reflected at the opposite wall 
and interact with the incoming upstream particles, forming a shock. We use 10000 x 
300 cells with a resolution Axi = Ax2 = O.SSc/wp, where oOp = 4:7rni2e'^ /me is the 
plasma frequency with the upstream electron density measured in the downstream 
frame at t = 0. The number of particles per cell is 3 x 3, the time step is 0.25 w"^, and 
the total simulation time is 4800 o;"^. 
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Figure 2. Shock structure at t — 2395 Wj^ : the charge density in the X1-X2 plane 
(a), its spatial average (b), the phase space diagram (c) and the B3 component of the 
magnetic field (d). 



The shock, which propagates along the positive xi direction, is formed after 
t ^ 350c<;~^ Figure |2] shows the important physical quantities at t = 2395 w^^. The 
typical filamentary structure of Weibel-mediated shocks ahead of the shock front can be 
seen in the charge density as well as in the magnetic field (figure |2^, d) and the density 
compression factor is ~ 3 (figure [2|d). The phase space diagram (figure [2]:) shows the 
thermalized downstream region on the left hand side and the shock transition region 
with escaped and reflected particles on the right hand side of the shock front (located 
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around 1000 c/ujp for the conditions of Fig. |2| 

At early times, the filamentary structure does not affect significantly the shock 
structure and its influence at later times is addressed in the following section. However, 
averaging over the transverse spatial component gives good qualitative agreement 
between the theoretical estimates and the simulation results throughout the entire shock 
propagation. 



4. Discussion 



The analysis of the density ratio associated with the shock front shows that this ratio 
can reach up to ^2/^12 = 3.2 ± 0.08. This illustrates that when the shock structure 
is generated self-consistently the shock density ratio can deviate from the theoretical 
value n2/ni2 = 3, which is derived from the jump conditions for a cold plasma [20] and 
a Maxwellian distribution in the downstream with a thermal spread A7 = 9.5 (leading 
to r2 = 1.525). We also observe a slight deviation from the shock velocity /3g2 = 0.49 
{^measured = 0.48). Siuce the impact on the density ratio is clearer, we will limit our 
detailed discussion to this quantity. In order to analyze the impact of the accelerated 
particles on the jump conditions we have measured the adiabatic constant directly from 
the kinetic information of the particles in the simulation data as well as analytically 
from the fittings to the data in figure [3^. For the analytical estimate we assume a 



particle distribution given by equation (11 ). Both methods provide essentially the same 
results. The adiabatic constant decreases logarithmically from initially r2 = 1.5258 to 
r2 = 1.5247 at the end of the simulation (figure [3b) which predicts a density change 



according to equation (A. 4) of 6n2/ni2 ~ 0.01 and does not explain by itself the density 
deviation which we observe in the simulations. We note that the changes in the adiabatic 
constant are very small and the fluctuations of the data points are almost on the same 
level as the total decrease in F. 
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Figure 3. (a) Electron distributions for simulation times tcup ~ 900 (blue), 4800 (red) 
with Maxwellian fit with A7 = 9.5 (dashed) and indication of a power-law with index 
a = —2.9 (solid black), (b) Evolution of the adiabatic index for electrons (blue) and 
positrons (red) and logarithmically decaying fit. 
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The particle distribution function in the downstream region is almost homogeneous 
along Xi and varies slowly, whereas the physics in the shock transition region is highly 
dynamic. In the following, and in order to calculate the pressure and charge densities 
along the shock propagation direction, the particle distribution ahead of the shock is 
treated as a single bulk stream, which might not be appropriate for large simulation 
times, but in the early stages (up to tujp ^ 2000), the fraction of escaped or reflected 
particles is low compared to the bulk. The pressure density profile along xi is used to 
define the integration range for the quantities ahead of the shock, the Lorentz factor 712 
and the upstream enthalpy /ii. The peak in the pressure is considered as the transition 
between upstream and downstream regions and the integration range is varied up to 
300 c/ojp. 




500 1000_H 1500 2000 

time [u ] 




Figure 4. (a) Specific entlialpy (blue) and average Lorentz factor (red). In light 
colors, the contribution of positrons (dashed) and electrons (solid) is shown. The 
integration region is 100 c/wp ahead of the shock front, (b) Normalized field energies 
^i/^kin with i — (red), Ei (black), E2 (blue) and currents ji (green), j2 (dashed). 
El is multiplied by 5. 

After the shock is formed, the Lorentz factor ahead of the shock deviates strongly 
from the initial value 7^2 = 20 (Fig. |4^), which leads to an increase of the density 
ratio according to Figure [TJ At the same time, the specific enthalpy has increased, 
which has a decreasing effect on the density ratio. Both quantities are oscillating in 
phase, where a high enthalpy appears together with a low bulk Lorentz factor and vice 
versa. The decrease of the average Lorentz factor in front of the shock stems from a 
mixing of different populations ahead of the shock: the incoming upstream flow with 
'j^2 = 20, which is decelerated by the fields at the shock front, the particle precursor, 
which consists mainly of escaping particles that have not been affected by the shock, and 
the reflected particles from the upstream region. Our simulations reveal that scattering 
of the flow impinging on the shock front in the self-consistent fields generated in the 
shock front leads to significant heating (in both the longitudinal direction and in the 
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transverse direction) at the expense of the free energy of the flow, thus contributing to 
the overall slowdown of the flow as it approaches the shock front. 



]^ (Time = 500 m^) \ , (Time = 1000 (1)^ ) 

-2 -10 12 -2 -10 12 




Figure 5. 2D plots: Longitudinal and transverse currents in the shock region at times 
tujp = 500, 1000. ID plots: Average fields \Ei\ (black), IB3I (green), \E2\ (gray); 712 
(red). 

Fig. |4]d shows that the electric field component \Ei \ grows while the Lorentz factor 
is decreased and reaches a first saturation point at ~ 60 w"^. The growing magnetic 
field converts energy from the longitudinal momentum pi to the transverse component 
P2, which causes charge separation between positrons and electrons. The associated 
current ji is responsible for the appearance of Ei, showing its peak value at the same 
time when the longitudinal field \Ei\ saturates after the linear stage. As the average 
value <Ei> is zero, statistical changes in the longitudinal electric field component must 
be responsible for the slowdown of the particles in the shock region. At tujp > 500 the 
magnetic field B3 and the transverse electric field E2 are increased and the positron 
and electron species start to oscillate in antiphase around a mean value. E2 is smaller, 
but close to Bs and follows the same trend, which is typical for Weibel-type instability 
generated filaments. The transverse fields are generated via instabilities of Weibel-type, 
which generate and amplify fluctuations in the longitudinal currents. In Fig. |5]the total 
currents ji = ji^e+ + ji,e- {i = 1,2) are plotted at tUp = 500, when the transverse field 
components start to grow and the oscillations in the species become strong, and at 
tUp = 1000, when the quantities in Fig. have reached a quasi-steady state. While ji 
is strong in the entire region of the particle precursor and very weak behind the shock, 
j2 exists only in a sharp region around the shock front, and coincides with the peak in 
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Table 1. Definition of the models shown in figure [6j The bullets indicate if 
the deviations from the adiabatic constant = r2 — 5T2 and the Lorentz factor 
7^2 — 712 + 5^12 are taken into account and lS.xtrans denotes the transition region. 
The bulk is a Maxwellian in all models. 



Figure [6] compares the average downstream density from the simulation with the 
different theoretical models listed in table [H It is clear that the simulation results 
differ from the ideal model (Ml - no changes in 712 and enthalpy). The inclusion 
of deviations from a Maxwellian distribution function of the downstream (M2) does 
not affect the density ratio significantly. On the other hand, the contribution of the 
decreasing upstream Lorentz factor is observed to have an important impact on the 
density ratio, but strongly depends on what is defined as the upstream region of the 
shock. The comparison of the results for different integration ranges shows that after an 
initial overshoot, the quasi steady state solution of the jump conditions for an integration 
range of lOOc/wp (M5/M6), matches the data best. This suggests that only the vicinity 
of the shock front within this range significantly affects the shock properties. 

3.4 

3.3 

CM 3.2 
c 

3.1 



3.0 



500 1000 _i 1500 2000 

time [m ] 

Figure 6. Downstream density from simulation data (solid black) and comparison 
with theoretical models according to table [ij Ml - dashed black, M2 - orange, M3 - 
dashed red, M4 - dashed blue, M5 - red, M6 - blue. 




If the contributions from the self-generated electromagnetic fields are considered. 
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the resulting density ratio is slightly decreased. In equation (|4]), the first term on the 
left-hand side and the pressure term on the right-hand side are the dominant terms, and 
of the same order {p2/'n2U2s ~ uufii ~ 20). For the magnetization to become important, 
let us assume a contribution of 10%, so that it has to exceed a = 0.05 as f3is ~ 1. The 
total magnetization in our simulations is ~ 0.05 after a quasi-steady state has been 
reached, which makes it necessary to be included in the discussion of unmagnetized 
shocks. The additional decrease of the density ratio due to this contribution is of the 
order of 0.1, which is calculated from the conservation equations ([T])-(|4]). 

5. Conclusions 

In conclusion, we have investigated the evolution of the shock properties when 
corrections from the usually considered fluid theory are taken into account due to the 
self-consistent evolution of the shock. We have shown that the shock jump conditions 
are affected by these corrections, in particular the density ratio. We found that 
the formation of a non-thermal tail in the particle distribution and the associated 
monotonous decrease of the downstream adiabatic constant, as well as the modifications 
of the upstream bulk speed directly in front of the shock, lead to an increase of the 
density ratio. The build-up of the upstream pressure and electromagnetic fields have a 
decreasing effect on the density ratio. 

Results from 2D particle-in-cell simulations confirm our theoretical predictions, 
showing a density ratio 7% larger than predicted from the standard jump conditions, 
for early propagation times. The evolution of the upstream Lorentz factor (which has 
been demonstrated to slow down when approaching the shock ^) is observed to be 
the main quantity responsible for such deviations. This analysis allowed us to define 
the spatial range that determines the shock transition region, which is observed to be 
lOOc/cjp, illustrating that the shock is mainly determined by the particles and fields 
within this range. Our results open the way for a more detailed understanding of the 
self-consistent evolution of the shock properties, where kinetic effects are taken into 
account, and demonstrate that a quantitative comparison between shock parameters 
and simulations/observations should take into account deviations from the standard 
jump conditions. 
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Appendix A. Changes in the jump conditions due to the Lorentz factor 



The full expression of the shock speed is given by 

A. = (A.1) 

MiV7i2 - 1 

To obtain the influence of the upstream Lorentz factor 712 = 712 — ^li2, "we do a Taylor 
expansion, yielding 



o _oO (r2-l)(7i°2-/^l) . 

Ps2 — Ps2 \ — 7~n — "^712- 



From the density ratio 



we obtain 



ni2 f3s2 7i2(r2 - 1) (712/^1 - 1) 



n^^nl /ii(l-2/ii7i°2 + (7i°2)') e 
ni2^ni2^(7?2)'(r2-l)(7iVi-l) 



(A.2) 



(A.3) 



(A.4) 



Appendix B. Analytical expressions for the energy and pressure densities 

For a distribution function consisting of a Maxwellian plus a power-law tail and an 



exponential cutoff, deflned by equation (11), the analytical expressions for the energy 



and pressure densities deflned in equations ^ and ( 10 ) are given by 
6 = 271 {CiA7(l + 2A7(1 + A7)) exp (-A7-^) 



+Q 



/ mm 



'~1cut 



a 



A7 



it" exp 



Icut 



p 



TT { 2CiA72(l + A7) exp(-A7-^) + C2 



T 2 



a 



7; 



2-a 



Icut 
'^cut 



In 



a 



a 



+^lcut exp 



''Icut 



A7 



A7c«tr 2 



a, 



''icut 



'~)cut 



CB.l) 



cut / L \ ^'~)cut J \ ^'Jcut , 

where F(n, stands for the incomplete Gamma function. The constant C2 is obtained 
from the simulations, whereas the normalization condition determines 



(27r 



,-1 



cut) 



[a - 1 



1-1 



(B.2) 



A7(l + A7) exp (-A7-1) 

with the exponential integral function Ej^z). 
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